Dynamical response of ultracold interacting fermion-boson mixtures: Fermion-hole vs 

polaron-hole excitations. 
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We analyze the dynamical response of a ultracold binary gas mixture in presence of strong boson- 
fermion couplings. Mapping the problem onto that of the optical response of a metal/semiconductor 
electronic degrees of freedom to electromagnetic perturbation we calculate the corresponding linear 
response susceptibility in the non-perturbative regimes of strong boson-fermion coupling using dia- 
grammatic resummation technique as well as quantum Monte Carlo simulations. Depending on the 
interaction strength and the details of the underlying bosonic spectra we find that energy absorb- 
tion/emission can occur due to excitation of individual fermions as well as fully dressed polarons. 
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I. INTRODUCTION 

Ultracold gas mixtures of bosons and fermions, which 
arise very naturally in the sympathetic cooling processes, 
are very interesting quantum systems with physical prop- 
erties very different from conventional quantum gasesJ^ 
In view of recent advances in the field of correlated ul- 
tracold gases it is very important to understand their 
dynamical properties, e. g. their response to an exter- 
nal dynamical scattering potential brought about by a 
change of trapping potential. If the bosonic subsystem 
is dominantly in the BEC phase the effective low-energy 
interaction with the fermionic subsystem is the coupling 
between the latter and the phonons (sound waves) of the 
former. From the mathematical point of view such a 
system is nothing but an electron-phonon coupled sys- 
tem best described by the Frohlich Hamiltoniani^ In case 
of low fermion concentrations it describes individual im- 
purities imbedded into a continuum of massless bosonic 
modes. Under such conditions the physics of the sys- 
tem is supposed to be very close to that of the classi- 
cal polaron, taking place in semiconductors with strong 
electron-phonon interaction4i"— 

There are, however, fundamental differences between 
the conventional (solid state) polarons and their BEC 
counterparts. The most obvious one is the different 
phonon spectrum of the bosonic subsystem as well as an 
explicit momentum dependence of the electron-phonon 
coupling)^ While these details do not alter the general 
picture of polaron static properties (there is still an effec- 
tive mass generation and self-trapping), they could possi- 
bly alter the dynamical response, which reveals such im- 
portant information as how the impurities interact with 
their surroundings In this paper we would like to con- 
sider them in full detail and in different geometries with 
the special emphasis on strong coupling results thereby 
closing the gaps in the existing literature. 

The paper is organized as follows. In the next sec- 
tion we formulate the problem and introduce all relevant 
quantities. Section Hill is devoted to the non-perturbative 



approach inspired by the classical random phase approx- 
imation (RPA) . We explain the details of the implemen- 
tation and discuss the special features pertinent to ul- 
tracold gas realizations. In Section IIVI the calculation 
of the dynamical response function is accomplished us- 
ing numerically exact quantum Monte Carlo simulation 
technique. Section fVl contains a discussion of results and 
offers several avenues of further progress. 



II. THE MODEL AND OBSERVABLES 

An effective low-energy Hamiltonian for a BEC- 
fermion mixture has the canonical Frohlich form^i^ 
which is written in terms of boson (described by the field 
operators &k) and fermion (denoted by Oq) degrees of 
freedom, 



q k 
q k^O 



(1) 



where the dispersion of the fermion is Eg = q'^ /2m, 

LUk = ck[i + my2]'/^ (2) 

is the dispersion of the phonon mode and the coupling is 
given by = A[(efc) V((efc)2 + 2)] V4 ^ith A = 51B V^Vq. 
giB is the effective interaction strength between impuri- 
ties and Bogoliubov excitations and can be adjusted by 
changing the particle density and/or the s-wave scatter- 
ing length of collision processes of the impurity with the 
bosonic medium. ^ denotes the healing length of the con- 
densate and is given by ^ = 1/ ^/STraesno where as b is 
the boson-boson s-wave scattering length and uq is the 
condensate density. 

One fundamental difference between the 'classical' 
semiconductor based electron-phonon coupled system 
and the one realized in ultracold mixtures is that the 
quantum gas system can be prepared in trapping poten- 
tials for fermions and bosons which might be of different 
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shape and dimensionahty. For that reason we shall later 
consider systems with different Eq and Wk- Changing the 
shape of the trapping potential for the impurity in space 
and time, for instance by acceleration with respect to the 
BEC, which rests in the laboratory reference frame one 
induces the rearrangement of particles. Another possible 
probing process is the Bragg spcctroscopyi^"— An appro- 
priate response function to that is the autocorrelation of 
the particle current density. This picture directly corre- 
sponds to the conventional polaron, where the response 
function is charge current autocorrelation, which is re- 
lated to the frequency-dependent or optical conductivity 
ct(w)J^ In order to estabhsh immediate parallels to the 
classical polaron problem we shall call the quantity in 
question optical conductivity as well. 

In general the optical conductivity is a complex quan- 
tity. While the imaginary part describes attenuation 
of the excitations, the real part is the actual emis- 
sion/absorption spectra. We shall concentrate on the 
latter. It can be shown that 

Re[cr(w)] = -Im[7r(w)/w] , (3) 

where 7r(cj) is referred to as force-force correlation func- 
tion, which is up to some frequency dependent prefactor 
proportional to the correlation function of the perturba- 
tion on the rhs of ((T])J^ The calculation is most conve- 
nient in the Matsubara representation, where we have 

g2 rf> 
^ (*'*^) ^ 2r ^2^/ / dre'"^ 

q,q' 

X {Trp (q, r) p (q', 0) $ (q, r) $ (q', 0)) (4) 
with the density operator 

P (q) = XI "k«k+q , (5) 

k 

and the generalized potential 

$(q) = K,(feq + 5t_J . (6) 

The subscript p, in the double sum indicates the spa- 
tial direction with respect to which the conductance is 
probed, that is p € {x,y,z}. After the computation of 
7r(icj) one has to perform an analytic continuation and 
then extract the optical conductivity using the prescrip- 
tion da]). 

Before proceeding with the actual calculation, a re- 
mark is in order. From the perspective of a solid-state 
physicist, the optical conductivity describes the exper- 
imental conditions quite well, that is probing a sample 
with optical or X-ray photons does not lead to a sub- 
stantial momentum transfer (Ap « 0). In case of the 
RF-spectroscopy in ultracold quantum gases, this is not 
necessarily the case. In order to use the abovementioned 
force-force correlation function, the absence of momen- 
tum transfer is essential. The Monte Carlo scheme dis- 
cussed later, however, is quite more general and allows 



for the consideration of an arbitrary momentum transfer. 
In this case, a more suitable observable is the real part of 
the fully momentum-dependent optical conductivity^^ 

R.e[cr^^ (q, w)] = Im[n^^ (q, w)] , (7) 

ui 

where, again, the most convenient form for the current- 
current correlation function n(q, w) is the Matsubara 
representation, 

n^. (q, ^) = -^ (Tril ^) 0)) ' (8) 
with the current densities 

.•(q)— ;^E(k + f)4«.+.. (9) 

In the analytic approach as well as in the Monte Carlo 
scheme we restrict ourselves to the case of one dimen- 
sional systems, hence skip the vector notation. Experi- 
mentally, this can be motivated by the use of quantum 
gases in reduced dimensions. 

If not explicitly stated otherwise we use the polaronic 
units. Distances are measured in ^, time in units of 
TO^^/fi and energies in units of /i^/(m^^). 



III. RPA RESULTS FOR THE RESPONSE 
FUNCTION 

Analytical approaches have been widely used in the in- 
vestigation of the effects of electron-phonon interactions 
in metals and semiconductors (for the latest review, con- 
sult e. g. Ref. [3). Very recently, due to formal similar- 
ities much effort has been made to adapt these methods 
to Bose-Fermi mixtures^ii^"— . Although the applicabil- 
ity and accuracy of these approximation schemes have 
in general to be questioned)^ their simplicity allows for 
an easy way to get first insights into the complex nature 
of electron-phonon interaction. Especially when it comes 
to the interpretation of numerical data, such analytical 
models proved to be of much use and importance. 

Wc begin our calculations with a simple perturba- 
tive treatment of the optical conductivity using the 
force-force correlation function. Mainly, wc are inter- 
ested in the effect of different energy-momentum rela- 
tions LUk and interaction potentials Vk- We are look- 
ing at the conventional Frohlich model with an LO 
phonon {ujkiVk) = (fio, A/ |/c|), the acoustic phonon 
model {ujk,Vk) — (|fc| , A/-\/|fcj), the small momenta 
BEC-polaron modelSi {ujk,Vk) = (c |fc| , A-\/|A^) and, of 
course, the full BEC-polaron model as presented in the 
previous section. 

The leading order contribution to the conductivity is 
found by replacing the expectation value with respect to 
the full Hamiltonian in Eq. ^ by an expectation value 
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for the non-interacting system. 
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7r(ia;) 



where 

is the Green's function of the phonon 
mode and x(g, iuj + iv) is the density-density correlation 
function (polarization loop) defined as 

X'°' (<?, r) = -{Tr piq, t) p{-q, 0))o . (11) 

The corresponding Feynman diagram is depicted in 
Fig. [TJ The calculation of the free phonon propagator 



X(°) (r)I?(«) (r) 



FIG. 1. Leading order contribution to the conductance. The 
wiggly line represents the free phonon propagator and the 
solid lines denotes the free Green functions of the impurities 
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FIG. 2. (Color online) Real part of the conductivity for the 
discussed models and several chemical potentials. Panel a) 
the LO model (Qq ~ 1), b) the acoustic phonon model, c) 
the small momenta BEG polaron model and d) the BEG po- 
laron model. The temperature is fixed to /3 = 10 and the 
chemical potential is varied: fj. = 0, 0.5, 1, 2 (black solid, red 
dot-dashed, blue dashed and orange long-dashed line). 



as well as the free polarization loop is straightforward, 
one obtains. 



1 



1 



— UJq IV + UJq 



fl2) 



and 



IV + Ek - ek+q 



where np denotes the Fermi-Dirac distribution function 
and ek — Ek — fi with the chemical potential fi. The re- 
maining Matsubara summation and the analytical con- 
tinuation iujn — > Lo+iO^ are straightforward. One obtains 
for the conductance, 

Re [cr (w)] = — — / -;— s [riB {(^q) - (w, + rsu)] 

r,s—± 

x\q\V^n^{eq)\ ^ (14) 

where ne is the Bose-Einstein distribution function. 
In Fig. [2] the conductivity is depicted for the above- 
mentioned models for several chemical potentials. The 
Frohlich model shows the well-known threshold be- 
haviour (i.e. (T w for u < flo). At finite tempera- 
ture, each considered model shows a divergent tendency 
for small frequencies. This is due to the combination of 
the coefficient l/oj^ and the factor consisting of Bose- 
Einstein distribution functions. In the limit a; — >■ one 
usually expects to obtain the DC conductivity in the case 
of the 'canonical' polaron system. This feature is also re- 
ferred to as the Drudc peak. It is the static response to 



a constant electric field and might indeed be rather large 
as compared to the rest of the spectrum. In the BEC 
case this Drude peak is directly related to the mobility 
of the impurity. In order to obtain the purely dynami- 
cal response it is sensible to subtract the w = result. 
Interestingly, for temperatures approaching zero where 
the TT-B can be approximated by ub (w) = —Q (— w) this 
tendency is strongly suppressed (see Fig. |3]). 
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FIG. 3. (Golor online) Real part of the conductivity for the 
discussed models and several chemical potentials. Same con- 
ventions as in Fig. [2] except for the approximation ne (i^) ~ 

The effect of increasing chemical potential is a shift 
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of the pcak-likc structure towards higher frequencies. In 
the LO model an additional feature can be observed: a 
double peak structure emerges (one peak at w w fio and 
another one ai uj Qq + fi). For non-trivial energy- 
momentum relations ujq this feature is washed out. This 
double-peak structure has a natural explanation. The 
optical response measures how easy it is to excite an 
electron-hole pair. The existence of the first threshold 
is trivial: if there is not enough energy to overcome fio 
than the response is suppressed. With growing a; one has 
to 'dig' deeper into the Fermi sea. The maximal energy 
for the electron-hole pair is then equal to the band depth, 
which in the present case is equal to fj,. 

In order to go beyond the perturbative analysis, sev- 
eral field theoretical approximation schemes are avail- 
able. We have chosen to use the random-phase approx- 
imation, which is good for systems with efficient screen- 
ing. Usually this happens only in systems at high den- 
sities. This assumption, of course, is questionable for a 
large fraction of experiments concerning ultracold quan- 
tum gases. For these systems (where Migdal's theorem 
is a priori not applicable), an approach based on vertex 
corrections seems to be more promising. Nonetheless, the 
present approximation is a good starting point for more 
sophisticated techniques. On the other hand, we shall see 
later that RPA captures many effects found by numer- 
ically exact Monte Carlo approach. Therefore here we 
focus on systems with sufficiently high particle densities. 

In our RPA we compute Eq. (j4|). 



quency summation. Then one obtains^i^ 



7r(ia;) 



(iuj)^V ^-^ 



(15) 

using the approximated phonon Green's function 
given by 



(16) 



The corresponding diagrammatic representation of the 
screened phonon propagator is depicted in Fig. |3) Due 

^VWW^ = "fx/x/N/x/v/^ -I- *\/\/\/\/\/^ >VWW^ 

T Q "^iX^^/^a 

FIG. 4. Screened electron-phonon interaction: Dyson-like 
equation of the phonon propagator in the random phase ap- 
proximation. 

to the more involved structure of I?^^^, the Matsub- 
ara summation with respect to iv is not simple anymore. 
Therefore, it is more suitable to use the Lehmann repre- 
sentation of the Green's functions to get rid of the fre- 



TT [iLo) 



x-^ oof dE 
(luj) V ^ J 

{2Im [pRPA.fl- (q^ E)] x(°' (g, lio + E) 



2Im 



2Im 



RPA 



V 



RPA 



{q,E)x'''Hl^E)Y} (17) 



where the superscript R denotes the retarded component 
of the corresponding Green's functions. Again, we are 
only interested in the real part of the conductivity, hence, 
the last expression in Eq. ([T7)) can be neglected. Because 
of the simple structure of the free polarization loop, the 
remaining integral can be performed yielding 



Re [cr (w)] 



s=± k,q 



Im[2?RPA'«(g,sefc+,g-.sefc 



X [ns (sEk+sq - SEk - W) - TIB {sek+sq " SCfe)] . (18) 

We restrict the calculation to the zero temperature case. 
In this situation the polarization diagram possesses a 
closed analytical expression. 



x(«.*-)--^^E'«"t-'f^^ll:|, (19) 



This allows for an efficient calculation of the screened 
phonon propagator. In the case of finite temperature 
and quadratic dispersion relation for the impurities, the 
evaluation of the polarization loop is numerically rather 
involved. Figs. [5] and |6] show the conductivity for several 
parameter constellations. 

As already seen in the perturbative calculation. Fig. [5] 
shows a pronounced shift of the secondary peak with 
changing chemical potential. The most prominent fea- 
tures, however, are seen for changing coupling strength, 
see Fig. [51 In the case of the LO phonon the positions of 
the maxima in the two-peak structure are consistent with 
the perturbative calculation and are located at w = J7o 
and = f^o + In the case with linear phonon dis- 
persion relation ~ for the acoustical phonons and a BEG 
system with linearized dispersion [panels b) and c) ] the 
position of the secondary peak is different. In fact its 
location is still very sensitive to changing but without 
any perceivable shift at different A. 

Another interesting feature is the pseudogap, seen as 
the first minimum in the spectrum, which develops in the 
non-BEC case. Its width grows with coupling strength 
and is directly connected to the energy stored in the fully 
developed polaron state which has to be destroyed prior 
to excitation of the impurity. In the BEG case, no pseu- 
dogap formation is observed. The reason for that is, of 
course, not the absence of the polaron state as seen in 
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interacting systems. In order to corroborate the RPA 
results and clarify the physical picture one has to em- 
ploy more advanced techniques. Quantum Monte Carlo 
simulation method is one of such powerful approaches 
which enable to access the dynamical properties of the 
system i^^i^^ Although in some cases it is subjected to the 
limitations such as finite size and the sign problem, this 
method has been successfully used in many polaron re- 
lated problems Its results have also been considered 
as benchmarks where exact solutions are not available. 

Our Monte Carlo simulation is based on a path inte- 
gral formulation of the dynamical correlation functions. 
We shall follow the theoretical treatment worked out in 
a previous work^S with an extension of it to the conju- 
gate momentum space which is convenient in our present 
study. 



FIG. 5. (Color online) Conductance the RPA approximation 
for fixed interaction strength A = 1 and varying chemical 
potential = 0.5, 1, 1.5, 2 (black solid, red dot-dashed, blue 
dashed and orange long-dashed line) . The order of the graphs 
is the same as in the figrues above. 
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FIG. 6. (Color online) Conductance in the RPA-like approx- 
imation for fixed chemical potential fi — 1 and varying in- 
teraction strength A = 0.25,0.5,0.75, 1 (black solid, red dot- 
dashed, blue dashed and orange long-dashed line). The order 
of the graphs is the same as in the figures above. 



the LO case. It seems that in the EEC case, due to the 
linearity of the bosonic dispersion, the polaron cloud is 
rather 'rigid' and the whole polaron can be excited to the 
higher lying fcrmion state (in the absorption case). On 
the contrary, the Drude peak tends to be enhanced. 



IV. MONTE CARLO SIMULATIONS 

As we have seen in the previous section, the opti- 
cal spectra change considerably in the case of strongly 



A. Formulation of path integral 

For our purposes it is very convenient to reduce the 
phonon field to the original set of harmonic oscillator 
operators as qu = 1/ ^2mpUJk{b^_^. + b^), and pk = 
i-\/mpajfc/2(&l^ ^bk), to replace the phonon creation and 
annihilation operators in Hamiltonian (1), with mp be- 
ing the effective oscillator mass. So that the Frohlicli 
Hamiltonian is rewritten as. 



Pi 



1 



q k 



(20) 



By using the standard Trotter's decoupling scheme, we 
represent the Boltzmann operator in a path integral form, 



VxTr exp <y — 

X{\xkm{xk{o] 



dT[he{T,x) + hph{T,x)] 



(21) 



where 

he{T,x) =^{Eq- fl) al {T)aq (t) 



EE 

q fc/0 



hp{T,x) = E ^ 



rUpUJk 



al+k{'^)aqij)qk{T) 



\ rup 


'dxkir)' 


I 2 





+ -jJUpU^kXt^ir) 



(22) 



Here r is imaginary time, \xk) is the eigenstate of the op- 
erator qk with eigenvalue Xk satisfying the eigen-equation 
qk\xk) = Xk\xk)- Then we define the time evolution op- 
erator Ux{t) along a path x as. 



(23) 
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On this path, the free energy (= F^) is evaluated by a 
trace over the fermionic part of the Bohzmann operator, 

e-PF^ = e-^"''^^''''(^)Tr[[/^(/3)] 

= e-^'''^^''''(^)det|I + U^(^)|, (24) 

where \JxiT) is the matrix representation of the time 
evolution operator Ux(t). The partition function (= Z) 
and total free energy is obtained by integrating out the 
bosonic field of x, 



Z 



(25) 



The expectation value of an operator O is given by 



(O) 



Vxe' 



where (O)^ is the average along a path x, 

TT[U4f3)0] 



Tr[C/,(/3)] 



(26) 



(27) 



In this notation, the path integral form for the single 
particle Matsubara Green's function is given by: 



G{q,T;q',T') = ^J Vxe'^'"'' G,{q, t; , 

G.,{q,T;q',T') = -{T,a,{T)a'^ {r')) (28) 

where ag(r) is the Heisenberg representation of Og, de- 
fined by 



(29) 



By solving the equation of motion for the involved 
operators^ we get the path-dependent Green's function 
(/? > T > r' > 0), 

G,((z,r;g',r') = - {U,(r)[l + U,(/3)]-i 

xU;i(r')}^,^,. (30) 

The optical conductivity and absorption spectrum are 
derived from the current-current correlation function ([SJ , 
which can be rewritten as 

n^i^g, T;q',T') = J 'Dxe~^^''Il^^^x{q,T;q' ,t'), 
^tj.i^.x{q,T;q',T') = -^{Trjl{q,T)ji,{q' ,t'))x 

' i:(''..+l)(*i+l 



kk' 



X [G,(fc', r' + /?; k, T)Gx{k + q, r; k' + g, r') 

+ G,(fc, I3;k + q, 0)G,{k' + q, /3; k' , 0)] . (31) 

In the last line, the time-dependent Bloch-Dc Dominicis 
theorem^S has been employed to decouple the many-body 
operators into a product of bilinear components. 



B. Numerical results on dynamical responses 

In our numerical calculation, the path integral is per- 
formed by the quantum Monte Carlo simulation method 
along with the matrix factorization and QDR decom- 
positions in quad precision to reduce the numerical er- 
rors. During the data acquisition, the dynamical correla- 
tion function is measured after about every 100 updates 
and we take totally 5000-20000 samples until the simu- 
lation converges. In addition, extra 100-500 samples are 
swept to thermalize the starting configuration to equi- 
librium regime before the acquisition. After evaluating 
the current-current correlation function via Monte Carlo 
simulation, we derive the absorption spectrum R^^(ci,uj) 
by inverting the integral equation 



(32) 



n^iy(g, t) = -/ dujRf^^{q,uj 



1 - e-^"' 



From Rf^^{q,uj), the optical conductivity is easily ob- 
tained as before by taking the following limit 



lim —Ri_,^{q,uj) 



(33) 



In order to solve the integral equation Eq. ([5^ for the 
two-particle correlation function, we have developed a 
renormalizing iterative fitting method. Some details of 
this algorithm are described in the Appendix. 
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FIG. 7. (Color online) Real part of optical conductivity of 
a ID conventional LO model with 25 polaronic states when 
A=0.5, 13=10 and Qq = 1. The three spectra in the main panel 
correspond to filling fractions: single particle (black solid), 5 
(red dot-dashed) and 11 (blue dashed) particles. The inset 
displays the same results in semi-logarithmic scale to expose 
the fine structure. 

In Fig. [3 we plot the real part of optical conductivity 
a of Frohlich model calculated from the quantum Monte 
Carlo simulation. Our simulation is performed on a ID 
lattice with 25 sites with open boundary conditions. The 
total number of fermions (= nf) is adjusted between 1 
and 25 by changing the chemical potential. Since we can 
only deal with models of finite size, here and after we con- 
fine our simulations in a momentum space of — tt < q < tt 
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(here measured in inverse lattice spacing), where the po- 
laronic states are uniformly distributed. In Fig. [71 we 
plot the spectra for three typical cases of 5 and 11 

particles. In addition to the Drude peak at around cj=0, 
the spectra show a phonon peak at cj=2. 5^3.0 coming 
from the single-phonon process. In order to visualize the 
fine structures in the spectra, we replot these data in the 
inset with a semi-logarithmic scale. Here, one can clearly 
identify that besides the Drude and phonon peaks, there 
is another island of weak excitations in the spectra. This 
island is due to the high-order electron-phonon scattering 
processes and is located at the high energy part. When 
we increase Uf from 1, we actually introduce more po- 
larons into the system. So that the many-body effect 
emerges gradually when the polarons begin to interact 
with each other. As a result, in the figure we observe 
that with Hf increase, the single-phonon peak is some- 
what broadened, and the multi-phonon island is shifted 
towards left with an enhanced amplitude. This tendency 
suggests that the polaron effect can be reinforced by in- 
creasing the density of fermonic impurities. 



ing temperature one can see the single-phonon peak gets 
closer to characteristic energy oi lu = Hq. 




FIG. 8. (Color online) Temperature dependence of real part 
of optical conductivity for a ID conventional Frohlich model 
with 11 impurities on a lattice with 25 sites: p—20 (black 
solid), /3=10 (red dot-dashed) and /3=5 (blue dashed). The 
coupling constant is fixed to A=0.5. The inset shows the spec- 
tra in semi- logarithmic scale. 

As illustrated in Fig. [S] from RPA calculation, the 
single-phonon peak is located around u;=1.0, that is at 
the frequency fig of the LO phonon. In Fig. [71 however, 
the single-phonon peak appears at a higher energy. This 
is a finite temperature effect. At T = 0, the Fermi sea 
has a sharp edge at the Fermi energy Ep. Therefore, 
in the single-phonon process, the creation of an electron 
just above Ep and a hole below Ep has a sharp thresh- 
old at Hq. At finite T the Fermi edge is washed out, but 
symmetrically both towards smaller and larger energies. 
On the contrary, the bosonic part is deformed in a dif- 
ferent fashion, which leads to a perceived energy shift of 
the threshold for the phonon excitation. To quantify that 
we have performed simulations at different temperatures, 
see Fig. [51 adjusting /? between 20 and 5. With dccreas- 
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FIG. 9. (Color online) Real part of optical conductivity of 
a ID conventional Frohlich model with an LO phonon with 
11 impurities on a lattice with 25 sites at /3=10. The three 
spectra in the main frame correspond to different coupling 
strength A: 0.5 (black solid), 0.7 (red dot-dashed) and 0.85 
(blue dashed). The inset displays the same results in semi- 
logarithmic scale. 

Next we investigate the dependence on the interaction 
strength, varying A between 0.25 and 0.85. The last value 
is of the order of the band width and thus we expect it 
to drive the system deep into the strong coupling regime. 
Fig. [9] shows the results for the LO phonon situation, 
again on a lattice with 25 sites and 11 impurities. All 
curves share the general feature of a peak around to ~ 2.3, 
which is just the usual threshold frequency given by the 
phonon frequency plus the level spacing due to the finite 
size. This compares well with the RPA calculation. The 
secondary peaks due to excitation of impurities which lie 
deeper in the Fermi sea are much less pronounced, but 
undergo a shift towards lower frequencies for growing A. 
This is much better seen in the inset of the figure. 

More fundamental differences can be observed for the 
Drude peak. Starting with A = 0.5 it is moving towards 
finite energies while the spectral weight at a; = almost 
completely vanishes. Thus a pseudogap opens up indi- 
cating that the system becomes 'insulating'. This effect 
is fully covered by the RPA calculation and reflects the 
polaron binding energy. 

Now we turn to the case of the BEG with the disper- 
sion relation ([2]). We first study its dependence on the 
density of fermonic impurities. In Fig. llOi the real part of 
optical conductivity of a ID 21-sitc system is depicted for 
filling fractions n/=:l/21, 5/21 and 11/21, respectively. 
Here the coupling constant is set to A=0.25, and the in- 
set displays the spectra in a semi-logarithmic scale so as 
to amplify the fine structures. From the black curve of 
n/ = l, one can infer that A=0.25 is already of intermedi- 
ate strength because there is a clear pseudo gap at uj=0. 
This is different from the RPA calculation. Thus we con- 
clude that the formation of the gap requires a presence 
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FIG. 10. (Color online) Real part of optical conductivity for a 
ID system involving fermions and bosonic excitations in BEC, 
when A=0.25, p—lO. The three spectra in the main frame cor- 
respond to different doping levels: single (black solid), 5 (red 
dot-dashed) and 11 (blue dashed) fermonic impurities. The 
inset displays the same result in a semi-logarithmic represen- 
tation. 



FIG. 11. (Color online) Real part of optical conductivity for 
a ID model with 11 fermionic impurities in a system with 21 
sites at 13— 10. The three spectra in the main frame correspond 
to different coupling strength A: 0.25 (black solid), 0.35 (red 
dot-dashed) and 0.45 (blue dashed). The inset shows the same 
result in semi-logarithmic representation. Blue arrows in the 
figure indicate the phonon shoulders. 



of an additional energy scale in the system, fig in the 
LO case or finite level spacing in the present situation. 
Due to numerical restrictions wc could not test out this 
hypothesis with sufficient precision though. 

With Uf increasing to 5, one finds that the gap width 
and depth change, as shown by the red dot-dashed 
curves. If we go on increasing to 11 as well as for larger 
fillings (not shown), as illustrated by the blue dashed 
curve, the gap eventually closes. Since the level spacing 
of an interacting system usually tends to decrease with 
growing number of particles this is consistent with the 
picture offered above. 

There is also an alternative explanation for the pseu- 
dogap formation though. It might be related to the k- 
dependence of phonon energy w/c = ck[l + {S,k)'^ /2]^^^ . As 
already known, the phonon modes of momenta k ~ i:2q-p 
are responsible for the gap opening, where qp is the Fermi 
momentum. When we raise up E'f, Qf and uj2qj, also in- 
crease. If u!2qp is very large, these phonon states would 
become prohibitively difficult to generate. So that the 
gap will not appear. Such evolution of the gap cannot 
be observed in the RPA calculation, which is done for 
the continuum model from the outset. However, as the 
state-of-the-art apparatus allows for generation of rather 
short optical lattices we believe that the pseudogap could 
be observed experimentally. 

The above mentioned difficulty to open a gap at high 
n f can be partially eliminated by increasing the coupling 
strength A. The attraction between fermions and bosons 
produces a negative coupling energy and tends to com- 
pensate the energy gain from phonon creation processes. 
In the ID system, this coupling would end up with a gap 
opening, provided that A is large enough to stabilize the 
phonon states of 2qp. In Fig. Ill[ wc examine this prop- 
erty from a view of optical conductivity for the BEC po- 



larons. We again consider 11 fermions immersed in a ID 
trap of BEC with 21 sites. The gap appearance is simi- 
lar to that in Fig.lH Here A=0.25 (black sofid), 0.35 (red 
dot-dashed) and 0.45 (blue dashed) correspond to weak, 
intermediate and strong couplings, respectively. How- 
ever, in Fig. 1111 the phonon peak is highly modified not 
only in its shape but also in its position. When A=1.5, 
the phonon peak even plunges into the Drude peak, giv- 
ing rise to a broad shoulder. This phonon shoulder is de- 
noted with an arrow in the main graph as well as in the 
inset. Such a strong dclocalization effect on the phonon 
peak is absent in the conventional Frohlich model. It can 
be attributed to the /c-dependence of coupling Vk in the 
BEC polaron model. The most important features of the 
spectra are concentrated around a few k values. Since 
Vk = A[(C/s)V((^fc)2 + 2)]i/4, ff A increases to a larger 
value A', we can find a smaller k to keep Vk invariant, 
i.e. Vk'{X') = Vfe(A), thus applying a 'discrete' version 
of rcnormalization transformation. That means wc can 
attain roughly the same coupling energy by enhancing A 
and reducing k together. Therefore, the phonon excita- 
tions gradually concentrate to the small k regime, and 
the low energy phonon modes become more favorable as 
A increases. This displacement of phonon peak leads us 
to conclude that for large coupling, those phonon modes 
of small momenta or long wavelengths play more impor- 
tant roles in the dynamical response of the system. 



V. DISCUSSION AND CONCLUSIONS 

We have analyzed the spectrum of the correlation func- 
tion of particle currents in a number of interacting mix- 
tures of bosons and fermions. We considered different 
kinds of couplings and dispersion relations for the con- 
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stituent subsystems, which model electrons in semicon- 
ductors coupled to LO phonons, acoustical phonons as 
well as fermionic impurities which are immersed into a 
BEC and which interact with Bogolyubov modes of the 
condensate. While for the solid state realizations the 
quantity we calculate is the optical response, in the BEC 
case the correlation function of currents gives a direct 
access to the Bragg spectra of impurities. 

While in the weak coupling case we recover all of the 
known physics, we find distinct effects in the strongly in- 
teracting case. Especially a pseudogap (with respect to 
particle-hole pair excitation) formation could be found 
using both the RPA approach as well as Monte Carlo 
simulations in ID. We speculate that this effect is due 
to the finite energy stored in the polaron state, which 
is released/absorbed during the fermion excitation pro- 
cess. Remarkably, this is pertinent only to solid state 
realizations with an LO phonon and BEC systems of a 
finite size. Large BEC systems does not show up such 
behaviour due to linearity of their bosonic spectra in the 
low-fc sector. Therefore in this regime the energy ab- 
sorbing/emitting process occurs in such systems without 
polaron destruction. 

Our approaches allow for an extension to a number of 
realistic experimental setups. We expect that the effects 
we predict could, for example, be investigated in binary 
systems which are similar to those used in Refs. [3l| or 
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Appendix A: Analytic continuation of two-particle 
Green's function 



In order to calculate the optical conductivity, we have 
to derive the absorption spectrum uj) from the in- 

tegral equation 



11^1^(1?, t) = - / dujRfj^t,{q,uj 



which is called analytic continuation. In our numeri- 
cal calculation, we use a renormalization iterative fit- 
ting method to do this. The original iterative fitting 
method is developed for the one-body Green's function 
of electrons^. It relies on the sum-rule of electronic spec- 
trum, which conserves the total spectral weight through 
the iteration process. However, for the two-body spec- 
tra like absorption spectrum R^^{q,ui), there is no such 
a rule that conserves the total spectral sum like that of 
electron, except for the following properties. 



Ri_,y{q,~u) = -R^i,{q,u). 



(A2) 
(A3) 



Obviously, the two-body spectrum is anti-symmetric and 
the total sum of spectral weight should be zero. 



dujR^y{q,Lo) = 0. 



(A4) 



Therefore, the standard iterative fitting method does not 
work here. 

Nonetheless, because the spectrum is anti-symmetric, 
we confine our calculation in the region < cj < 00. We 
also introduce a modified spectral function 



Rp.u{q,i^) = - 



np.(q,/3) 



coth ( ^— 



(A5) 



Because of the anti-symmetry of Rfi^{q, w), on substitut- 
ing R^^{q,ijj) in Eq. (jAip . we can rewrite the integral 
equation into 



1 fP^ 



da;i?^y(g, w) cosh ( — 



X cosh ( ) np,.(9,/3). 



(A6) 



Here R^i,{q, w) is just a rcnormalizcd form of the original 
spectral function, but it is easy to see this new function 
is positive for w > 0, and it satisfies a sum rule 



du>Rf^„{q,uj) 



1, 



(A7) 



which allows us to solve the integral equation of Eq. (jA6[) 
with the iterative fitting algorithm in the regime < 
Lu < 00. Once R^u{q,Lo) is obtained, the original spectral 



1 — e^*^" 



, (Al) function Ri^^{q,uj) can be reproduced from Eq. (jA5 
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